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We present a simple and efficient method to simulate three-dimensional, complex-shaped, in- 
teracting bodies. The particle shape is represented by Minkowski operators. A time-continuous 
interaction between these bodies is derived using simple concepts of computational geometry. The 
model (particles -I- interactions) is efficient, accurate and easy to implement, and it complies with 
the conservation laws of physics. 3D simulations of hopper flow shows that the non-convexity of the 
particles strongly effects the jamming on granular flow. 

PACS numbers: 02.70.Ns 45.70.-n 45.40.-f 47.ll.Mn 



Molecular Dynamics, MD is the art of modeling com- 
plex systems as a collection of particles interacting each 
other. MD in three-dimensions using arbitrary particle 
shape is a fundamental problem in several disciplines: 
Drug molecules often act as a key in a lock formed by a 
protein cavity, so that they can be designed using MD 
simulations Liquid crystals consisting of complex- 
shaped molecules exhibit transition to a nematic phase, 
which can be investigated from particle-based simula- 
tions of complex shaped objects 0; The large scale mod- 
eling of geological materials in foundations, landslides 
and fault zones requires a constitutive equation which 
can be constructed using MD-like models |3|; Comput- 
ing the motion of rigid and articulated bodies can lead 
to new advances in robotics, automation Q and virtual 
reality applications 

The most typical approach for these applications is 
to solve the dynamics of interacting rigid bodies, where 
their real shapes are approximated by polyhedra [1, 0, 
01. The most difficult aspect for the simulations is to 
model contact interactions. Contact force methods have 
been proposed for two-dimensional (2D) models using 
polygons^ Q. However, the extension of this method 
to three-dimensional (3D) simulations has proven to be 
extremely difficult. In the simple case of convex poly- 
gons, the force is calculated as a function of their overlap 
area However, the assumption that clastic force is a 
function of their_overlap leads to a non-conservative elas- 
tic interaction 8] . An alternative approach is to assume 
that the potential elastic energy is a function of the over- 
lap. Then forces and torques are derived from this poten- 
tial 's'l . Both approaches have still not been extended to 
3D, because the calculation of the overlap between two 
polyhedra is computationally very expensive. This is the 
main reason why most of the commercial codes for partic- 
ulate systems are still based on simulations with spheres, 
or clumps of spheres representing complex shaped objects 

i- 

An alternative solution for the 2D simulations of com- 



plex shaped particles has been proposed recently 
The method introduces the concept of spheropolygons, 
which is the object resulting from dilating a polygon by 
an sphere. The method not only guarantees energy bal- 
ance but also proves to be much more efficient that previ- 
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ous models to represent complex particle shape 
In 3D models, the dilation of a polyhedra by a sphere has 
a precise mathematical meaning using the Minkowski op- 
erator. In our knowledge, Liebling and Pournin were the 
first to introduce the Minkowski operators in particle- 
based simulations 
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In order to calculate the in- 
teractions, they assumed a single contact point between 
the particles [l3|. These approach, however, leads to 
forces discontinuous in time and numerical artifacts such 
as abrupt creation of mechanical energy. An alterna- 
tive approach is proposed by Pournin by calculating the 
overlap area between the particles [l5^ . This approach in 
practice is not feasible due to the high complexity of the 
boundary of the sphcropolyhedra. 




FIG. 1: The spherotetrahedron (right) is obtained by sweep- 
ing an sphere into a tetrahedron (left). 



In this Letter we present a solution to this problem by 
using the concept of spheropolytopes. They are gener- 
ated from the Minkowski addition of a polytope by an 
sphere, which is nothing more than the object result- 
ing from sweeping a sphere around a polytope. A poly- 
tope is a generic mathematical concepts that can refer 
to polygons, polyhedra or polygonal curves in 3D. The 
polytope is regarded as a collection of features in the 
three-dimensional Euclidian space: vertices, edges and 
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faces. The interaction between two spheropolytopes is 
calculated as a function of the distance between their 
features. The molecular dynamics is implemented in a 
simple, efficient and elegant algorithm, which complies 
with conservation laws of physics. We believe that this 
model will lead to a wide range of applications of molec- 
ular dynamics, as complex particle shape and realistic 
interactions can be captured in a unified framework us- 
ing well established concepts of molecular dynamics and 
computational geometry. 

For the representation of arbitrary shaped particles we 
introduce the mathematical concept of Minkowski sum. 
Given two sets of points P and Q in an Euclidean space, 
their Minkowski sum is given hy P(BQ = {x + y\ x& 
y ^ Q}- This operation is geometrically equivalent to 
the sweeping of one set around the profile of the other 
without changing the relative orientation (Fig[T]). Exam- 
ples of Minkowski sums are the spheropolygons (sphere 
® polygons) [lot , spherocyllinder (sphere © line seg- 
ment) [l3|, the spherosimplex (sphere © simplex) 11411 
and the spheropolyhedron (sphere ® polyhedron) [15l |. 
All these objects can be enclosed in the generic shape 
of spheropolytopes, which consists of a set of vertices 
Vi, edges Ei, polygonal faces Fi and the radius r of the 
sweeping sphere wich will be called sphero-radius. The 
three examples of spheropolytopes we will consider in this 
letter is the rice, the tetra and the yermis in the Fig. [2l 




FIG. 2: Spheropolytopes generated as sphere line segment 
{rice), sphere © tetrahedron ( tetra) and sphere © polyline 
{yermis). 



imal of the distance between each vertex of one edge and 
the other edge. 

For the calculation of the distance between a vertex 
Vi and and face Fj we consider the plane Ilj containing 
the face. First we project the point on the plane. If the 
projection of the point in the plane lies inside the polygon 
face then the distance is calculated between these two 
points. If the projection lies outside the polygon then the 
distance is calculated to the closest point in the polygon 
boundary. 

These formulas of distance are used to calculate the 
force between two spheropolytopes. the force Fij on the 
i-spheropolytope by the j-spheropolytope is taken as a su- 
perposition of the interaction between each pair of edges 
F{Ei,Ej) and each pair of vertex-face F{Vi,Fj) for the 
spheropolytope pair. 



Ei.Ei 



(1) 

The force F{Gi,Gj) associated to the two features (edge- 
edge or vertex-face) is assumed to depend on the over- 
lapping length 5 between them 



5{Gi, Gj) — Ri + Rj — d{Gi, Gj), 



(2) 



with d{Gi,Gj) the distance between the features of 
the spheropolytopes Ri the spheroradius of the i-th 
spheropolytope. The point of contact between the two 
features is calculated by taken the spheres of radius Ri 
and Rj centered in the closest points Xi and Xj, and 
finding the intersection between the line connecting these 
two points and the line connecting the two intersection 
points of the spheres. This contact point results as 



For the calculation of the contact force between 
spheropolytopes, we require expressions for the distance 
between their features. Given two features Gi and Gj of 
the spheropolytopes i and j, their distance is defined as 
d{Gi, Gj) = \ \Xi — XjW, where Xi and Xj is the closest 
points belonging to either feature. We start with the for- 
mula of the distance between a vertex Vi and an edge Ej . 
Let's consider the line ii containing the edge. First we 
calculate the closest point of this line to the vertex, if 
the point lies on the edge this is the closest point of the 
edge to the vertex. Otherwise we take the minimal of the 
distance from the vertex to both endings of the edge. 

Next we consider the distance between two edges Ei 
and Ej. Let's the lines £i and £j containing the edges, 
we find the two points on these lines whose distance is 
minimal. If both points belong to the edges, they define 
the minimal distance between the edges. Otherwise the 
distance between the two edges is calculated as the min- 



R{Gi, Gj) = X 



Rjf — R^ + (P{Gi, Gj) 



(3) 



where fi 



Xi-Xi 



2d{G„Gj) 

. , ^ ^ . I . From the point of application of the 

\\Xj-Xi\\ 

contact forces we get the torque on the i-spheropolytope 
by the j-spheropolytope: 

T,, = E r(£;.,£;,)-f E ^(^^'^.'O+E ^(^j-'^')- (4) 

Ei,Ej Vi,Fj Vj,Fi 

where 



Ti{G,,Gj) = {R{G„G,) ~ r'^) X F{G„ ) (5) 

with fi the center of mass of the i-spheropolyhedron. 

Since the formulas of distance are continuous functions 
on the degrees of freedom of the spheropolytopes, the to- 
tal force is continuous too. This avoids the problems of 
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discontinuity in time of the forces in previous models [i| . 
Different forces can be included in this model: for ex- 
ample, a force derived from a potential function of the 
distance leads to a conservative systems; forces depend- 
ing on the relative velocities at the contact points leads to 
dissipative granular materials; Forces depending on the 
history of relative velocity at the contacts represent fric- 
tional granular systems; more sophisticated forces can be 
used to simulate biomolecules. The electrostatic interac- 
tion between the molecules can be modeled by allowing 
the force depend on the closest points between the fea- 
tures. 



FIG. 3; Spherotetrahedra collider for checking the conserva- 
tion laws. 
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TABLE I: Percentage error in the numerical simulations cal- 
culated from the mechanical energy (E), angular momentum 
(L) and linear momentum(p) refereed to their initial values. 
The simulation time is 20 s, the mass of the particles is 1. kg 
and the stiffness is 10000 N/m 



time step is smaller. The discretization error is larger 
in the yermis-particles, because they produce much more 
collisions than the other particles. 





The first numerical experiment presented in this pa- 
per is the sphcropolytopes collider (SPC). The simulation 
consists of colliding opposing beams of sphcropolytopes. 
We consider the three shapes shown in Fig. [5] . The 
SPC is build with the intention to produce enough colli- 
sions to test the basic conservative laws of physics: the 
conservation of energy, linear and angular momentum. 
The experimental setup is showed in Fig. [31 two rows of 
spherotetrahedra are set in a collision course with ran- 
dom orientation and angular velocity. The contact force 
between features is calculated as 

F{G,,Gj)^^kn6{G,,Gj)n, (6) 

where fc„ is the elastic constant. Once all the forces are 
calculated we integrate Newton's second law using the 
Verlet algorithm for the translation coordinates. The 
Euler equations form angular momentum is integrated 
using and the Fincham Leap Frog algorithm, based on 
the quaternion formalism, for the orientation coordinates 
The Table n shows the percentage error of the three 
conservation laws (energy, linear and angular momen- 
tum) for two different time steps. The consistency of our 
numerical method is verify as the error decreases as the 




FIG. 4: Hopper flow simulations at initial (left) and flnal 
(right) stages for rice (above), tetra (center) and yermis (be- 
low). The parameters are the same than in the SPC experi- 
ment. The coefficient of viscosity 7 is equal to 0.2 s"^ 



The immediate extension of the model is to include 
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visco-elastic forces: 



P{Gi, Gj) — —knS{Gi, Gj] 



(7) 



where 7 is the viscocity force and Vc is the relative 
vefocity of the particles at the contact. This contact 
force offers an interesting application of this model: the 
study of the effect of particle shape on the jamming 
phenomenon of granular flow. The flow may happen 
when particles are discharged through a small opening, 
but particles may became jammed when the opening is 
smaller than a critical value. Modeling of gravity flow 
has been done using circular or spherical particles 
but the effect of shape on flow has not been thought- 
fully investigated. In particular, non-convex particles is 
expected to jam more easily than convex, or circular par- 
ticles. 

Granular flow with convex and non-convex particles 
is presented using the same three particles geometries 
shown in Fig. [21 The simplicity of our model allow us to 
represent the hopper and the container as just another 
example of spheropolytopes, see Fig. IH Contrary to 
previous findings [8| our convex shaped particles do not 
become clogged in the hopper. This is because we have 
not introduced an static frictional force yet. However, 
as an striking result, the non-convex particles get stuck 
without static friction, as can be seen in Fig. [5l 
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FIG. 5: Number of particles exiting the hopper in the setup 
shown in Fig [l] The simulation consider 90 particles. Note 
that in the yermis case the flow stops at 12 s due to clogging. 



Modeling interacting particles using spheropolytopes 
has several advantages with respect to other existing 
particle-based models: i) The possibility to model non- 
convex particles (in our case yermis, hoppers and con- 
tainers) ; ii) a realistic representation of the surface curva- 
ture of particles; iii) guaranteed compliance with physical 
laws; iv) numerical consistency, guaranteed by the conti- 
nuity in the proposed contact law. v) efficiency, given by 



a simple model for the contact law relying only on dis- 
tance calculations. This is radically different from pre- 
vious approaches where the contact forces are calculated 
in base of overlaps 0, 

The most interesting aspect of this model is to pro- 
vide a general framework for generic particle shape and 
contact interactions. Spheropolytopes is a very general 
shape which can be uses to represent biomolecules, poly- 
mers, rocks, meteors, etc. For the modeling of geologi- 
cal materials, particles with random shapes and tunable 
roundness can be generated by applying Minkowski oper- 
ators on Voronoi diagrams |12l |. Comminution processes 
can be model by solving the continuum stress equation 
of a given spheropolyhedra and originate from it fracture 
planes and hence secondary spheropolyhedra. 

This work is supported by the Australian Research 
Council (project number DP0772409 ) and the UQ Early 
Carrier Research Grant. 



[1 
[2: 

[3: 
[4: 
[s: 

[6] 
[7] 

[8 

[9: 
[10: 

[11 

[12 
[13 

[14 
[15 
[16 
[17 



Electronic address: ls.galinotorres@uq.edu.au Also at 
Physics Department, Grupo de Simulacion de Sistemas 
Fisicos. Universidad Nacional de Colombia. 
Electronic address: ]fernando@esscc. uq.edu. au| 

H. Carlson and A. Mccammon, Molecular Pharmacology 
57, 213 (2000). 

G. Pelzl, S. Diele, and W. Weissflog, Advanced Materials 
11, 707 (1999). 

F. Alonso-Marroquin, S. Luding, H. Herrmann, and 

I. Vardoulakis, Phys. Rev. E 51, 051304 (2005). 

B. Mirtich, ACM Transactions on Graphics (TOG) 15, 
177 (1998). 

D. Ruspini and O. Khatib, Journal of Robotic Systems 
18, 769 (2001). 

D. Baraff, Algorithmica 10, 292 (1993). 

S. Hasegawa and M. Sato, Computer Graphics Forum 

23, 529 (2004). 

T. Poeschel and T. Schwager, Computational Granular 
Dynamics (Springer, Berlin, 2004). 

M. Lu and G. R. McDoweU, Granular Matter 9, 69 

(2007) . 

F. Alonso-Marroquin, Europhysics Letters 83, 14001 

(2008) . 

F. Alonso-Marroquin and Y. Wang (2008), 
arXiv:0804.0474. 

S. A. Galindo- Torres and F. Alonso-Marroquin (2008), 
submitted to Phys. Rev. E, arXiv:0811.2858vl. 
L. Pournin, M. Weber, M. Tsukahara, J.- A. Ferrez, 
M. Ramaioli, and T. M. Liebling, Granular Matter 7, 
119 (2005). 

L. Pournin and T. Liebling, in Powders & Grains 2005 
(Balkema, Leiden, 2005), pp. 1375-1478. 
L. Pournin, Ph.D. thesis, Ecole Polytechnique Federale 
de Lausanne (2005). 

Y. Wang, S. Abe, S. Latham, and P. Mora, Pure and 
Applied Geophysics 163, 1769 (2006). 
K. To, P. Lai, and H. Pak, Physical Review Letters 86, 
71 (2001). 



